/*******************************************************************************

	Authors: Federico Curci (federico.curci@cunef.edu)
	         Sebastien Fontenay (sebastien.fontenay@uclouvain.be)
			 Federico Masera (f.masera@unsw.edu.au)

	Version: 1.0
	Last update: Nov. 2018

*******************************************************************************/

cap program drop discretiz
program define discretiz
version 11.0
syntax varname [, range(string) interval(string) endogenous(string) depvar(string) second exogenous(string) interact(string) xt(string) vce(string) print save graph(string)]

quietly {

// Error messages

if "`range'"=="" {
	noisily display in red "[, range(min/max)] is missing"
	exit
}

if "`interval'"=="" {
	noisily display in red "[, interval(min(step)max)] is missing"
	exit
}

if "`endogenous'"=="" {
	noisily display in red "You must specify an endogenous variable with [, endogenous(varname)]"
	exit
}

if ("`second'"!="" & "`depvar'"=="") {
	noisily display in red "You must specify a dependent variable with [, depvar(varname)]"
	exit
}

if ("`second'"=="" & "`depvar'"!="") {
	noisily display in red "You must specify the option [, second] to estimate the second stage regression"
	exit
}

// Define macros and matrices to store results

local max=substr("`range'", strpos("`range'", "/")+1, .)
local j=.
if "`second'"=="" {
	matrix results=J(1,5,.)
}
else {
	matrix results=J(1,4,.)
}

// Discretize variable and estimate regression

forvalues i=`range' {
	forvalues k=`interval' {
		local j=`i'+`k'
		if (`j'<=(`max')) {
			local lb=subinstr("`i'", ".", "",  .)
			local ub=subinstr("`j'", ".", "",  .)
			tempvar treat_`lb'`ub'
			gen `treat_`lb'`ub''=(`varlist'>=`i' & `varlist'<=`j')
			if "`interact'"!="" {
				replace `treat_`lb'`ub''=`treat_`lb'`ub''*`interact'
			}
			** FIRST STAGE **
			if "`second'"=="" {
				if "`xt'"!="" {
					xtreg `endogenous' `treat_`lb'`ub'' `exogenous', `xt' vce(`vce')
				}
				else {
					reg `endogenous' `treat_`lb'`ub'' `exogenous', vce(`vce')
				}
				test _b[`treat_`lb'`ub'']=0
				matrix temp=J(1,5,.)
				matrix temp[1,1]=`i'
				matrix temp[1,2]=`j'
				matrix temp[1,3]=r(F)
				matrix temp[1,4]=_b[`treat_`lb'`ub'']
				matrix temp[1,5]=_se[`treat_`lb'`ub'']
				matrix results=results\temp
			}
			** SECOND STAGE **
			if "`second'"=="second" {
				if "`xt'"!="" {
					xtivreg `depvar' (`endogenous'=`treat_`lb'`ub'') `exogenous', `xt' vce(`vce')
				}
				else {
					ivregress 2sls `depvar' (`endogenous'=`treat_`lb'`ub'') `exogenous', vce(`vce')
				}
				matrix temp=J(1,4,.)
				matrix temp[1,1]=`i'
				matrix temp[1,2]=`j'
				matrix temp[1,3]=_b[`endogenous']
				matrix temp[1,4]=_se[`endogenous']
				matrix results=results\temp
			}
		}
	}
}
matrix results=results[2...,....]
if "`second'"=="" {
	mat colnames results=lb ub fstat beta se
}
else {
	mat colnames results=lb ub beta se
}
cap matrix drop temp

// Save results

if ("`save'"!="" | "`graph'"!="") {
	preserve
	clear
	svmat results, names(col)
	gen cilow=beta-1.96*se
	gen cihigh=beta+1.96*se
	if "`save'"!="" {
		save discretiz_results, replace
	}

// Graph F-statistics

	if "`graph'"=="fstat" {
		gsort -fstat
		local lbval=lb[1]
		local ubval=ub[1]
		scatter lb ub [w=fstat], title("F-statistics") ytitle("Lower bound") xtitle("Upper bound") msymbol(circle_hollow) mcolor(green) name(fstat, replace) yline(`lbval') xline(`ubval')
	}

// Graph coefficients

	if "`graph'"=="coef" {
		gen nsign=(cilow<=0 & cihigh>=0)
		sum nsign, meanonly
		local nsignmean="`r(mean)'"
		gen possign=(beta>0 & nsign==0)
		sum possign, meanonly
		local possignmean="`r(mean)'"
		gen negsign=(beta<0 & nsign==0)
		sum negsign, meanonly
		local negsignmean="`r(mean)'"
		gen absbeta=abs(beta)
		if "`second'"=="" {
			local title="First stage coefficients"
		}
		else {
			local title="Second stage coefficients"
		}
		if ("`negsinmean'"!="0" & "`possignmean'"=="0" & "`nsignmean'"=="0") {
			graph twoway (scatter lb ub [w=absbeta] if negsign==1, legend(label(1 "- and sign. coeff") on) msymbol(circle_hollow) mcolor(blue))  , title("`title'") ytitle("Lower bound") xtitle("Upper bound") name(coef, replace)
		}
		if ("`negsinmean'"=="0" & "`possignmean'"!="0" & "`nsignmean'"=="0") {
			graph twoway (scatter lb ub [w=absbeta] if possign==1, legend(label(1 "+ and sign. coeff") on) msymbol(circle_hollow) mcolor(red))   , title("`title'") ytitle("Lower bound") xtitle("Upper bound") name(coef, replace)
		}
		if ("`negsinmean'"=="0" & "`possignmean'"=="0" & "`nsignmean'"!="0") {
			graph twoway (scatter lb ub [w=absbeta] if nsign==1  , legend(label(1 "n.s. coeff") on)        msymbol(circle_hollow) mcolor(yellow)), title("`title'") ytitle("Lower bound") xtitle("Upper bound") name(coef, replace)
		}
		if ("`negsinmean'"!="0" & "`possignmean'"!="0" & "`nsignmean'"=="0") {
			graph twoway (scatter lb ub [w=absbeta] if negsign==1, legend(label(1 "- and sign. coeff")) msymbol(circle_hollow) mcolor(blue))  ///
						 (scatter lb ub [w=absbeta] if possign==1, legend(label(2 "+ and sign. coeff")) msymbol(circle_hollow) mcolor(red))   , title("`title'") ytitle("Lower bound") xtitle("Upper bound") name(coef, replace)
		}
		if ("`negsinmean'"!="0" & "`possignmean'"=="0" & "`nsignmean'"!="0") {
			graph twoway (scatter lb ub [w=absbeta] if negsign==1, legend(label(1 "- and sign. coeff")) msymbol(circle_hollow) mcolor(blue))  ///
						 (scatter lb ub [w=absbeta] if nsign==1  , legend(label(2 "n.s. coeff"))        msymbol(circle_hollow) mcolor(yellow)), title("`title'") ytitle("Lower bound") xtitle("Upper bound") name(coef, replace)
		}
		if ("`negsinmean'"=="0" & "`possignmean'"!="0" & "`nsignmean'"!="0") {
			graph twoway (scatter lb ub [w=absbeta] if possign==1, legend(label(1 "+ and sign. coeff")) msymbol(circle_hollow) mcolor(red))   ///
						 (scatter lb ub [w=absbeta] if nsign==1  , legend(label(2 "n.s. coeff"))        msymbol(circle_hollow) mcolor(yellow)), title("`title'") ytitle("Lower bound") xtitle("Upper bound") name(coef, replace)
		}
		if ("`negsinmean'"!="0" & "`possignmean'"!="0" & "`nsignmean'"!="0") {
			graph twoway (scatter lb ub [w=absbeta] if negsign==1, legend(label(1 "- and sign. coeff")) msymbol(circle_hollow) mcolor(blue))  ///
						 (scatter lb ub [w=absbeta] if possign==1, legend(label(2 "+ and sign. coeff")) msymbol(circle_hollow) mcolor(red))   ///
						 (scatter lb ub [w=absbeta] if nsign==1  , legend(label(3 "n.s. coeff"))        msymbol(circle_hollow) mcolor(yellow)), title("`title'") ytitle("Lower bound") xtitle("Upper bound") name(coef, replace)
		}
	}
	restore	
}

// Display results matrix

if "`print'"!="" {
	if "`second'"=="" {
		mata : st_matrix("results", sort(st_matrix("results"), -3))
		mat colnames results=lb ub fstat beta se
	}
	else {
		mata : st_matrix("results", sort(st_matrix("results"), 4))
		mat colnames results=lb ub beta se
	}
	noisily matrix list results
}

}
end
